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We introduce a mean-field and perturbative approach, based on clusters, to describe the ground 
state of fermionic strongly-correlated systems. In cluster mean-field, the ground state wavefunction 
is written as a simple tensor product over optimized cluster states. The optimization of the single¬ 
particle basis where the cluster mean-field is expressed is crucial in order to obtain high-quality 
results. The mean-field nature of the ansatz allows us to formulate a perturbative approach to 
account for inter-cluster correlations; other traditional many-body strategies can be easily devised in 
terms of the cluster states. We present benchmark calculations on the half-filled ID and (square) 2D 
Hubbard model, as well as the lightly-doped regime in 2D, using cluster mean-field and second-order 
perturbation theory. Our results indicate that, with sufficiently large clusters or to second-order in 
perturbation theory, a cluster-based approach can provide an accurate description of the Hubbard 
model in the considered regimes. Several avenues to improve upon the results presented in this work 
are discussed. 
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I. INTRODUCTION 

Despite some substantial recent progress, an accurate 
and efficient description of the ground and excited states 
of low-dimensional strongly-correlated fermionic systems 
represents an open problem in condensed matter physics 
and quantum chemistry. A common feature in strongly- 
correlated systems is the collective behavior displayed by 
fermions in low-lying states. Accordingly, approaches 
based on composite particles have been proposed for 
treating these systems. One notorious example is the 
resonating valence bond as a ground state candidate for 
high-Tc superconductors, as suggested by Anderson [l|. 

In this work, we use composite many-fermion cluster 
states to describe the ground state of strongly-correlated 
systems. Here, a cluster is simply a subset of all avail¬ 
able single-fermion states that we group (generally using 
a criterion of proximity in real space). We presume that 
an accurate zero-th order description of the ground state 
of the full system can be prepared as a product of cluster 
states, each being many fermion in character. Two key 
aspects in obtaining an accurate description are: 1) the 
many-fermion state in each cluster is determined in the 
presence of other clusters, and 2) the single-particle basis 
used to determine the grouping into clusters is fully op¬ 
timized. This optimization could in principle break the 
real space localization criterion but in practice it gen¬ 
erally does not. The resulting cluster mean-field (cMF) 
state is, by construction, guaranteed to provide a varia¬ 
tional estimate of the ground state energy that is lower 
than Hartree-Fock (HF), i.e., the standard mean-field of 
single-particles. The optimization provides not only the 
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optimal cMF state, but also a renormalized Hamiltonian 
expressed in term of cluster states. Traditional many- 
body approaches can then be used, on this renormalized 
Hamiltonian, to account for the missing inter-cluster cor¬ 
relations. 

Our work is inspired by McWeeny ii , who first con¬ 
sidered the properties of wavefunctions written as a ten¬ 
sor product of the state of several subsystems (or groups), 
which are mutually orthogonal. McWeeny realized that 
the density matrix of cluster product states can be easily 
expressed in terms of the density matrices of the individ¬ 
ual clusters. In related work, Isaev, Ortiz, and Dukelsky 
U considered, in their hierarchical mean-field (HMF) ap¬ 
proach, a similar ansatz to ours for the 2D Heisenberg 
Hamiltonian. We note, nevertheless, that the generaliza¬ 
tion to fermionic systems of the HMF approach used in 
Ref. 0 is not straightforward if the full Fock space within 
each cluster is treated. An attempt was performed in 
Ref. 0 to split the Fock space in each cluster according to 
its number parity: states with even (odd) parity where 
treated as bosonic (fermionic) degrees of freedom. An 
ansatz for the ground state was constructed in Ref. 0 as 
a tensor product of the bosonic and fermionic degrees 
of freedom; this decoupling, however, may not be justi¬ 
fied in all cases and can potentially result in unphysical 
states. 

Our approach differs from that used in Ref. 0, aside 
from the application to fermionic systems, in not re¬ 
quiring the individual clusters to share the same ground 
state. That is, the ground state of each cluster is opti¬ 
mized independently allowing for (translational and spin) 
symmetry-broken solutions. In addition, we here con¬ 
sider Rayleigh-Schrodinger perturbation theory (RS-PT) 
@ to second-order as a means to obtain a correlated ap¬ 
proach defined in terms of clusters. 

A closely related cluster product approach was also 
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proposed by Li Q . In that work, the ground state of each 
cluster was, nonetheless, not optimized in the presence 
of other clusters. The author did, on the other hand, go 
beyond perturbation theory into a coupled-cluster ansatz 
(so-called block-correlated coupled cluster (BCCC)) as a 
way to correlate the cluster product state. The BCCC 
approach has been used with high success in quantum 
chemistry to describe strongly-correlated molecular sys¬ 
tems using either a complete active-space @ or gener¬ 
alized valence-bond m reference states. 

A cluster product state is naturally connected with 
the tensor network (TN) techniques that have been gain- 
ing pop ularity for treating strongly-correlated systems 
[T2,[l2. In essence, the cluster product state is the sim¬ 
plest possible TN, a simple scalar (the bond or ancillary 
index dimension is set to 1), although the elements defin¬ 
ing the network are chosen as cluster states rather than 
single-particle degrees of freedom as often done. The 
consequence of using a scalar product is that the clusters 
become disentangled; more general TNs such as the ma¬ 
trix product states (MBS) used within the density matrix 
renormalization group algorithm (DMRG) intro¬ 

duce entanglement in the ansatz and can provide highly 
accurate solutions for strongly-correlated systems. The 
optimization of TN states beyond the simple MBS is, 
however, non-trivial and remains an area of active re¬ 
search p^ . 

Yet other wavefunction ansatze that are related to the 
cluster product states are the correlator product states 
(CBS) [l^ or entangled plaquette states (EBS) 

Here, a variational ansatz is expressed in terms of entan¬ 
gled cluster products, as opposed to a simple uncorre¬ 
lated product. The price to pay is that the evaluation of 
matrix elements becomes more cumbersome, and it often 
has to be carried out by stochastic means (within a vari¬ 
ational Monte Carlo framework). We note that Ref. [T^ 
proposed a non-stochastic algorithm for optimizing CBS. 

At this point, we want to clarify why we have decided 
to use simple cluster product states even when more pow¬ 
erful ansatze are already available (such as more general 
TNs or CBS). In our perspective, the power of a clus¬ 
ter mean-field approach has not been fully realized. In 
particular, symmetry breaking and orbital optimization 
can partially account for inter-cluster correlations (when 
expressed in terms of the original set of fermion states). 
Moreover, the fact that the cluster mean-field state con¬ 
stitutes the ground state of a mean-field Hamiltonian of 
which the full set of eigenstates can be easily constructed 
has often been overlooked. This allows us to formulate 
a perturbative strategy to introduce the missing inter¬ 
cluster correlations. Yet more powerful many-body ap¬ 
proaches (such as coupled-cluster theory) can be easily 
introduced as done by Li in the BCCC method. 

Our objective with this work is thus two-fold. First, we 
present the cMF formalism as well as provide details of 
the RS-BT formulation we use. (We refer to the second- 
order perturbative result as cBT2 in the remainder of this 
paper.) We describe in some detail the strategy used to 


optimizate the one-electron basis in which the cMF state 
is expressed. Our second objective is to apply these tech¬ 
niques to the simplest paradigm of st rong ly-correlated 
fermionic systems: the Hubbard model [Q in one- (ID) 
and two-dimensions (2D) in a square lattice. The ID case 
is exactly solvable [^, while the 2D model has been ex¬ 
tensively studied. We refer the reader to Refs. and 
for a survey of numerical methods applied to the 2D Hub¬ 
bard model. We compare our results with second-order 
unrestricted Mpller-Blesset (UMB2) (2^ and unrestricted 
coupled-cluster with singles and doubles (UCCSD) [^ . 
as well as with perturbative triples (UCCSD(T)), cal¬ 
culations. MB2 constitutes second-order RS-BT based 
on a HF wavefunction (using canonical orbital and or¬ 
bital energies), while coupled-cluster constitutes a non- 
perturbative approach that involves an infinite-order re¬ 
summation of diagrams. Our results show that cMF 
(cFT2) significantly improves upon HF (MF2) and can 
provide an accurate description of the ground state of the 
Hubbard model. 

The remainder of this work is organized as follows. In 
Sec. in we present the formalism behind cMF and cFT2. 
Section uni provides some practical computational de¬ 
tails regarding the calculations presented in this work. 
In Sec. m we present the results of cMF and cFT2 cal¬ 
culations for the half-filled ID and 2D Hubbard model, 
as well as for the lightly-doped 2D regime. A brief dis¬ 
cussion following the results is presented in Sec. El along 
with some ideas as to how to improve the calculations 
here presented. Lastly, Sec. ED is devoted to some gen¬ 
eral conclusions. In Appendix we show higher order 
perturbation results in a small lattice, while in Appendix 
m we discuss the applicability of cMF to weakly corre¬ 
lated systems. 


II. FORMALISM 

A. Hubbard model 

In this work, we focus our attention on the Hubbard 
model in one- and two-dimensions (in a square lattice). 
The Hubbard model [T^ describes a collection of elec¬ 
trons in a lattice (of finite size L) interacting through 
the Hamiltonian 

H = -tY^ (C> Cj.„ + h-c.) + D ^ ni,t (1) 

i 

where the notation (ij) implies interaction only among 
nearest-neighbors. Here, cj ^ {ci,a) creates (annihilates) 
an electron with spin a on site i of the lattice and 
rii^a = cl^Ci^cr. The Hamiltonian contains one-electron 
hopping terms and an on-site repulsion (U > 0) term. 
The hopping amplitude t is used to set the energy scale. 
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B. Cluster Mean-Field 

Consider a set of single-fermion states {\k), k = 
1,..., M}, where M is the size of the basis, that satisfies 
the appropriate set of boundary conditions for the sys¬ 
tem. These single-fermion states may be different than 
the ones used to define the problem; in the case of the 
Hubbard model, they may be obtained by a rotation of 
the lattice (on-site) basis states: 

\k)=al\-), ( 2 ) 

(3) 

We assume that a\ and ak satisfy standard anti¬ 
commutation rules (implying orthonormality of {|fc)}). 
A basis for many-fermion states can be constructed 
from properly antisymmetrized products of such single¬ 
fermion states. 

Let the single-fermion states be grouped, according to 
some criterion (such as proximity in real space), into clus¬ 
ters of size ■ ■ -Im where n is the number of such clus¬ 

ters. The Fock space in each cluster can be constructed 
using the single-fermion states that define it. Due to the 
orthogonality of single-fermion states in different clus¬ 
ters, the Fock space of the full system is simply given by 
the tensor product of the Fock spaces of all clusters. 

A second-quantized formulation in terms of cluster 
product states can also be established. Let A| ^ (A/^c) 
create (annihilate) the J-th many-fermion state in cluster 
c. This J-th state is a linear combination of many-fermion 
basis states {\fJ.)c} (possibly mixing states with differ¬ 
ent number of fermions) constructed as antisymmetrized 
products of the single-fermion states in the cluster. We 
formally write 

\I)o = A\j-),, ( 4 ) 

where |—)c is the vacuum state in cluster c. (We empha¬ 
size here that |—)c does not correspond to the state with 
no fermions in the cluster, but is simply a useful abstract 
construct.) 

An arbitrary state |'I') in the Fock space of the entire 
system can be formed as 

1^) = ' |^)n, (5) 

I J Z 

where c/i;j 2 ;...;Zn are linear coefficients. Here, the sum 
over / spans the full Fock space in cluster 1, and so on. 
Each state in the expansion above constitutes a cluster 
product state. Formally, each cluster product state is 
built as 

\I),\J)2 - ■ - AlJ-), ( 6 ) 

where |—) is an abstract vacuum state for the full system. 

In this work, we consider a cluster product (mean-field) 
state as a variational ansatz for the ground state wave- 
function. That is, the ansatz |<i)o) for the ground state is 


given by 

|<l>o) = |0)i|0)2---|0)„, (7) 

where we have indicated that the ground state (hence 
the 0 label) in each cluster is used to build the product 
state. The optimal cMF state is obtained by a variational 
minimization scheme, as outlined in Sec. EDI 

Having defined a ground state cluster product configu¬ 
ration, excited configurations can also be considered. We 
write them as 

|d>/.)=|0)i---|/),---|0)„, (8) 

|$«;J,)=|0)l---|/).---|J)r--|0)n, (9) 

\^n-,Jr.Kk) = |0)i • ■ • |/). • • • I J), ■ • ■ \K)k ■ ■ ■ |0)„, (10) 

for singly-, doubly-, and triply-excited clusters. A full 
configuration expansion can be written in terms of ex¬ 
cited configurations as 

I /#o 

+ EEE + ■ ■ ■ (11) 

i<j //o J/O 

This provides exact eigenstates for the full system. 

Before proceeding further, let us comment on the na¬ 
ture of the cluster product states considered in this work. 
We indicated above that the ground state of each cluster 
is expressed as a linear combination of the many-fermion 
basis states in it. That is, 

|0)c = E<.cIm)c, (12) 

n 

with dp ^ = c(/r|0)c, where ^ is a compound index of oc¬ 
cupation numbers in the subset of single-fermion states 
of the cluster. The expansion over states {|^)c} can be 
restricted, i.e., some of the {d^ coefficients may be set 
to 0. We note that if a given cluster state is expanded in 
terms of even- and odd-number parity (here referring to 
an even or odd number of fermions) states, commutation 
rules between the operators A| ^ and Aj^c' are not sim¬ 
ple, complicating the evaluation of matrix elements as 
hinted below. All the calculations included in this work 
restrict the expansion of cluster states to a given and 
(or, equivalently, n and TOs) sector within the cluster, 
but include the full Hilbert subspace with those quantum 
numbers. This was done in order for the cluster product 
state |<i)o) to be an eigenfunction of and N^. 

Consider the determinantal expansion of the full sys¬ 
tem. The exact ground state wavefunction is expressed 
as 

E ■ ■), (13) 

including all possible many-fermion states {|/i)c} within 
the cluster. In cMF, the coefficients in the expansion 
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above are not independent, but are parametrized accord¬ 
ing to 


- i dg 2 dg^a ' 


(14) 


This parametrization per mits us to put cMF in the con¬ 
text of TN states [HI, and CPS In cMF the 

coefficient of each determinant is parametrized as a scalar 
product of cluster states (with compound indices). That 
is, the cluster product state has intra-cluster correlations, 
but lacks inter-cluster ones. This is in contrast to a TN 
state, where ancillary or bond indices include explicit en¬ 
tanglement in the ansatz. While in cMF the compound 
indices v, etc. refer to different orbital subspaces, CPS 
use overlapping indices as a means of introducing entan¬ 
glement in the ansatz. The optimization of TN states 
and CPS is, nonetheless, more involved than that of cMF 
states. 

A cMF state constitutes a generalization of a single 
Slater determinant, and thus HF can be written as a 
cMF where the orbitals are grouped into clusters, as we 
next describe. It should be stressed that this is only pos¬ 
sible in the optimized single-particle basis (or a unitary 
rotation of it that does not mix particles with holes). The 
HF state is recovered in different ways: as a product of 
two clusters (one fully occupied and one fully empty), or 
as a product of M clusters (the holes being occupied and 
the particles being empty), or arbitrary related construc¬ 
tions. 

The model also contains other wavefunction ansatze 
commonly used in quantum chemistry such as the an- 
tisymmetrized product of strongly orthogonal geminals 
(APSG) [Hl - [27f . Here, each cluster would contain two- 
electrons in a subspace of orbitals that define each gem- 
inal. It also encompasses the multi-configuration self- 
consistent-field (MC-SCF) [28 


. [Hi m odel as well as the 
30l. 311 variant of it. The 


complete-active-space (CAS) 
latter can be considered as a three-cluster state: the core 
is fully occupied, the virtual set of orbitals is fully empty, 
and the state in the active space is expressed as an op¬ 
timized linear combination of all possible many-electron 
basis states in the appropriate Hilbert space. 


Here, for instance, 

= X! + {Ai\p)al ap) 

pGi,r^j 

X! (pq\v\rs)alalasar + (17) 

prGi,qs^j 

Hijki=-^ ^ {pq\v\rs) alalttsar + - (18) 

rGk,sGl 

Given that fermion operators {a|,,ap} act on specific 
clusters, the matrix elements can be evaluated straight¬ 
forwardly if all the cluster states have a well defined 
number parity (though care has to be taken to respect 
fermionic anti-commutation rules). For instance, 

ifpG2 at|$o) =±|0)iat|0)2|0)3--- , (19) 

where the sign depends on the number parity of |0)i. 
(The action of aj, (op) on a specific cluster can be easily 
expressed in the occupation number basis {|/r)c} within 
each cluster.) If |0)c is of mixed number parity, the eval¬ 
uation becomes more cumbersome. For instance, 

lip £2 a],|$o) =+ |0)^ oj,|0)2|0)3--- 

- |0)rat|0)2|0)3--- , (20) 

where |0)(^ denotes the even-number parity projection 
out of |0)i, i.e., |0)^ = |-k)(-f|0)i. 

We close this section by noting that, if the ground state 
in each cluster preserves number parity, then expectation 
values of single-fermion operators (such as c(0|a],|0)c, for 
p £ c) vanish. This further implies that all three- and 
four-cluster interactions vanish in (<l>g | H | $g). If the num¬ 
ber of fermion states within each cluster is fixed, then 
the expectation value (<I>g|i7|$g) can be fully expressed 
in terms of the one- and two-particle reduced density ma¬ 
trices within each cluster, as first noted by McWeeny Q. 
(Note that this implies that the cost of evaluating the 
energy of a cMF state scales as 0{n^), where n is the 
number of clusters used.) 


C. Matrix Elements 

We now turn to the evaluation of matrix elements over 
cMF states. The (two-body) fermionic Hamiltonian is 
expressed in second-quantized form (in the basis that de¬ 
fines the clusters) as 

^ ^ '^{pq\v\rs)ala\as ar, (15) 

pr pqrs 

where (p|t|r) are one-body and {pq\v\rs) are two-body 
integrals. The Hamiltonian can be expressed as a sum of 
single, two, three, and four-cluster interactions: 

H = E -f 'y ) Hjj -|- y ) y ) (i6) 

i i<j i<j<k i<j<k<l 


D. cMF Optimization 

In this section we discuss how the cMF state is opti¬ 
mized, that is, how the ground state |0)c in each cluster 
c is determined. We use a diagonalization strategy akin 
to that used in HF or multi-configuration self-consistent- 
field (MG-SGF) methods. 

The optimal set of coefficients {dgcl Eq- US) can 
be found by minimization of the energy subject to the 
constraint that the state |0)c remains normalized: 

^($o|i7|$o)-eo.c<e = 0. (21) 

where eg^c is introduced as a Lagrange multiplier. The 
above equation can be cast as an eigenvalue equation 
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that, at the same time, defines a zero-th order Hamilto¬ 
nian within the cluster, z.e., 

= ( 22 ) 

The ground state of the cluster Hamiltonian is obtained 
as its lowest energy eigenvector. (Note that this also 
gives eo,c the physical meaning of the energy in cluster 
c.) The cluster Hamiltonian can be found trivially. As an 
example, if all cluster ground states {|0)c} have a fixed 
number of fermions, is given by 

^ {p\i\r)al Or + ^ X! (wl'Cks)ap a, as ar 

pr^c pqrs^c 

+ Psqi{pq\v\rs) - (pg|i)|sr)). (23) 

pr£c c' qs 

Here, = c'(0|aj as|0)c' is the one-particle density ma¬ 
trix in cluster c'. Because the cluster Hamiltonian 
depends on the ground state density matrices of other 
clusters, the equations must be solved self-consistently. 
This represents a generalization of the MC-SCF method, 
where there are several active subspaces (the clusters) 
each with its own multi-configurational expansion. 


E. Orbital Optimization 

As discussed in the introduction, in order to realize the 
full capability of cMF states it is necessary to include 
the optimization of the single-fermion basis in which the 
grouping into clusters is defined, which we refer to as 
an orbital optimization. Otherwise, a cMF state may 
yield an energy that is even above HF, despite having 
significantly more flexibility in the ansatz. We describe 
in this section how this is accomplished in our work. We 
note that the orbital optimization in cMF states is akin to 
the same process performed in traditional MC-SCF (and 
CAS, by extension) calculations in quantum chemistry. 

Given the single-particle basis {|fc)} in which the cMF 
state is constructed, we aim to rotate this to a new basis 
{|fc)} in order to lower the energy. We relate the two 
basis by a unitary transformation (parametrized as the 
exponential of an anti-Hermitian operator), 

= exp(k) flfe exp(—k), (24) 

k ='^{Kpqalaq-h.c.). (25) 

p<q 

In particular, we define an energy functional 

E[k] = ($o|exp(-k)ff exp(k)|$o), (26) 

where the (complex) elements {Kpq} serve as variational 
parameters. 


With the optimized cMF state |<i)o) at hand, we can 
compute the gradient with respect to orbital rotations at 
K = 0 (i.e., the gradient evaluated at zero-rotation) as 


r - dE 

~ Ok* 


= - $ 


K—0 


[H, 


aq Op] 


4>r 


Similarly, the Hessian can be constrncted as 

A B 


with 


An 


H = 


d'^E 


B* 


dn*. Oku 

-i(*„ 


Bpq^rs — 


[[H,a\ap\,al a^] 
+ q,P ^ r,s, 
d^E 






1 


+ 5<‘I>0 


Ojt Ot,! , Qq dn 


[[H,u,q u,pj , u,g 

+ q,p -H- s, r. 


$r 


(27) 


(28) 


(29) 


(30) 


The gradient (and Hessian) can be used to find a direc¬ 
tion of energy lowering with respect to orbital rotations. 
Several comments are in order: 


• The energy and Hessian can be evaluated following 
the strategy described in Sec. Ill Cl for the evalna- 
tion of matrix elements. It can be shown that the 
cost of building the full gradient and Hessian scales 
as 0(n^) and O(n^), respectively, where n is the 
number of clusters in the system. 

• Once a direction of energy lowering is found (defin¬ 
ing a non-zero k), we perform a finite rotation (an) 
of the Hamiltonian integrals. The energy functional 
is re-parametrized with respect to orbital rotations 
in terms of the new single-particle basis. 

• The above strategy (re-parametrizing the energy 
functional at each step) is necessary as the evalna- 
tion of the orbital gradient (and Hessian) is not as 
simple when k ^ 0. This also prevents us from us¬ 
ing a quasi-Newton strategy to perform the orbital 
optimization. 

• Care has to be taken of handling linear dependen¬ 
cies between the orbital rotations and the coeffi¬ 
cients in each cluster expansion. For instance, if 
a cluster state |0)c is expanded in terms of a full 
Hilbert (or Fock) subspace, then orbital rotations 
within the snbset of orbitals that define the cluster 
c do not lower the energy. 

In this work, as we use a full Hilbert subspace to de¬ 
scribe each cluster, the orbital gradient and Hessian for 
intra-cluster rotations are not considered as degrees of 
freedom. 
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F. Perturbation Theory 


In standard RS-PT we aim to solve for the eigen¬ 
states of the Hamiltonian H given the simpler Hamilto¬ 
nian for which all eigenstates are known. In RS-PT, 
the second-order correction to the ground state energy is 
evaluated as 




\Vo,\^ ^ 

So ~ Sfi’ 


(31) 


where V = H — and = ($o|)^Im)- Here, /r la¬ 
bels the eigenstates of Hq and are the corresponding 
eigenvalues. 

The cMF state, as outlined in Sec. ttmi provides a nat¬ 
ural zero-th order Hamiltonian of which all eigenstates 
can be easily constructed. This is expressed as a direct 
sum of the zero-th order Hamiltonians of each cluster 


ijo ^H° + H° + ... (32) 

The eigenstates of such Hamiltonian are given by 

J )2 ■ • • \Z)r,,=en,J2,--,Zn\I)l\J)2 ' ' ' \Z)n (33) 
Sll\J 2 \— \Zn = £1,7 + £2,7 + ' ' ' + £rt,Z- (34) 


As described in Sec. Ill Cl H has up to four-cluster in¬ 
teractions, while is, by construction, single-cluster 
in character. If a full Hilbert subspace in each clus¬ 
ter is used, matrix elements between the ground state 
|$o) and singly-excited (cluster) configurations vanish 
due to a generalized-Brillouin condition. Therefore, only 
two, three, and four-cluster interactions contribute to the 
second-order energy. The evaluation of the corresponding 
matrix elements can be carried out in a similar fashion 
as the evaluation of (‘I’oli^l^’o)- Naturally, computing 
the four-cluster interactions is the most expensive step 
in evaluating the second-order energy, with a computa¬ 
tional scaling of 0{n*) in the number of clusters. 

As described in Sec. Ill HI in this work we have chosen 
to use cluster ground states which preserve the number of 
t- and 4,-electrons. In that case, several two, three, and 
four-cluster interaction channels can be identified in the 
Hamiltonian, as summarized in Tab. HI The cMF ground 
state |<i)o) interacts with excited cluster configurations 
following these channels. 

At this point, we clarify that in this work the zero-th 
order cluster Hamiltonian is used, without any modifi¬ 
cation, to generate the full Fock space within the clus¬ 
ter. It is possible to tweak the definition of the non¬ 
interacting Hamiltonian {e.g., by adding a level-shift) in 
specific Hilbert space subsectors in order to improve the 
convergence properties of the perturbation series. 


III. COMPUTATIONAL DETAILS 

The cMF and cPT2 calculations presented in this work 
were carried out with a locally prepared code. Most of 


the results use an unrestricted cMF (U-cMF) formalism, 
where f-orbitals are allowed to have a different spatial 
distribution than 4,-ones. Some of the results in ID lat¬ 
tices use a restricted (R-cMF) formalism, where the spa¬ 
tial distribution is required to be the same. Real orbitals 
are used in both cases. In all the calculations we use the 
same number of 'f- and 4,-orbitals in each cluster, which 
we denote as I and refer to as the size of the cluster. The 
number of f- and 4^-electrons in each cluster was held 
fixed (thus preserving n and ms within each cluster). [32j 
Although not enforced from the outset, R-cMF calcu¬ 
lations resulted in spin singlet eigenstates within each 
cluster. 

The full relevant TOs sector of Hilbert space within each 
cluster was used in constructing the cluster ground state 
|0)c. For small cluster sizes, the ground state in each clus¬ 
ter was found by a standard diagonalization of the local 
cluster Hamiltonian. For larg er cluster sizes, a Lanczos 
[ 3 ^ or a Jacobi-Davidson [3413^ algorithm was used to 
solve for the ground state. 

The orbital optimization was carried out using a 
pseudo Newton-Raphson approach. After optimizing the 
cluster mean-field state, a Newton step was taken in the 
direction of energy lowering (using the orbital gradient 
and Hessian). A finite rotation provided a new single¬ 
particle basis in which the cluster mean-field was reopti¬ 
mized. These two steps were alternated until convergence 
was achieved in both the cMF state and the orbitals. 
This is akin to the most common methods of optimizing 
MC-SCF wavefunctions in quantum chemistry (Sbl - fs^ . 
We note that a full Newton-Raphson approach (with the 
mean-field and the orbital optimization carried out con¬ 
comitantly) should be preferred [s^, but we have not 
used it in this work, [d^ A globally convergent algorithm 
was used to guarantee that the variational cMF energy is 
reduced in each orbital optimization step. As described 
in detail below, for 2D lattices several local minima can 
be found in the orbital optimization process. We have 
not attempted to use an algorithm to locate the global 
minimum. 

In U-cPT2 calculations, all relevant cluster states were 
used in computing the second-order energy for small clus¬ 
ter sizes {1 = 2 and 1 = 3). For I = 4 and I = 5, the 
four-cluster contributions were computed using only 16 
states in each Hilbert subspace of a cluster, while two- 
and three-tile contributions used all available states. In 
I = 6 calculations we truncated the number of states 
in each Hilbert subspace in three- (four)-tile interactions 
to 64 (16), while no truncation was done in computing 
two-tile interactions. An energy-based criterion for the 
cluster states was used to carry out the truncation. We 
should point out that the second-order energy appears to 
be converged in all cases with respect to the number of 
states included. 

UCCSD and UCCSD(T) calculations were carried out 
with a locally modified version of the MRCC code 
[4lj,|4^. Exact solutions to the ID Hubbard lattice were 
obtained by solving the Lieb-Wu equations. 
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TABLE I. Two-, three-, and four-cluster interaction channels when the ground state in each cluster has well defined n-f^ and nj,. 


# clusters type^ 

2 one-electron CT 

2 two-electron, opp spin, CT 

2 two-electron, same spin, CT 

2 two-cluster spin flip 

2 two-cluster dispersion 

3 two-electron, opp spin, CT 

3 two-electron, same spin, CT 

3 single-cluster spin flip 

3 one-electron CT -I- dispersion 

4 two-electron scattering 


sample interaction^^ restrictions 


a' 

q,a 

^s,a' 


pqs € i, r € j 



^P,CT 

4 - 

aO-s- 

a ^r,(T 

pq sr £ j 



^p,cr 


Cf'SjCT C 

r,f7 

pq € i, sr £ j 



^p,CT 

i- 


a CLr,iT 

ps £i, qr £ j 



^p.,(T 


^s,a' 

^r,(7 

pr £ i, qs £ j 



^p.,(T 

4 - 

qO-s- 

cr (^r,iT 

pq£i, s£ j, 

r e 

k 

^p,(T 

al,q 


V,t7 

pq£i, s£ j, 

r e 

k 

^P,(T 

4 - 

a 0'S- 

a ^r,(T 

ps£i, q£ j, 

r e 

k 


0/ 

q,a- 

^s,a' 

f^r,CT 

pr £i, q£ j, 

S £ 

k 


0/ 

q,C7 

^S,CT^ 

^r,(7 

p£i, q£j, s 

£ k, r £ 1 


CT denotes charge transfer. 

^ Only two-body interactions are shown. Here, a generic two-fermion interaction takes the form aj, ^ “r.a, with a' = cr or —a. 


IV. RESULTS 

In this section we present results of cMF and cPT2 cal¬ 
culations on the ID and 2D Hubbard models. We start by 
providing an illustrative example in Sec. IIV A1 where we 
get into some practical details regarding the optimization 
of cMF states and the way in which other results are pre¬ 
sented. In Sec. IIV 51 we consider the ID half-filled case, 
for which exact solutions are available. We then proceed 
to study the 2D half-filled case in Sec. IIV Cl and finally 
consider the lightly-doped 2D case in Sec. IIVDI Our 2D 
results are compared to highly accurate numerical esti¬ 
mates from Refs, and [i^ 


A. Illustrative example 

In this section we discuss some practical aspects re¬ 
garding the optimization of cMF states. In this way, we 
hope that the results presented in subsequent sections 
will become more transparent to the reader. We con¬ 
sider a 12-site Hubbard ID periodic lattice at half-filling 
and U/t = 4. For U-cMF calculations, we typically start 
from an unrestricted HF (UHF) solution; we take the 
resulting orbitals and perform a Boys localization [dj. 
Figure [T] displays, in the top-left scheme, (localized) oc¬ 
cupied and virtual spin-orbitals mostly tied to two sites, 
which are then used to define a cluster of 2-)" and 2-1 
orbitals in which a single electron of each spin is placed. 
The optimized state in the cluster is expressed as a linear 
combination of the 4 possible resulting configurations. 

After setting the initial orbitals and the corresponding 
tiling scheme {i.e., defining how orbitals and electrons 
are grouped into clusters), the cluster mean-field state 
is optimized by a self-consistent diagonalization of the 
appropriate cluster Hamiltonians. The orbital gradient 
(and possibly the Hessian) is then evaluated which deter¬ 
mines how orbitals in different clusters should be mixed 
in order to lower the energy. A new orbital basis is de¬ 


fined by, e.g., a steepest-descent step, and the cluster 
mean-field step is reoptimized in such basis. The pro¬ 
cess is repeated until convergence is reached in both the 
orbitals and the mean-field state. The top-right scheme 
in Fig. [T] shows the converged orbitals that define a sin¬ 
gle cluster in the calculation. In particular, the orbitals 
displayed are the natural orbitals (those that diagonalize 
the t" and 4-- one-particle reduced density matrix) mostly 
tied to the original two sites. The orbitals defining the 
cluster remain well localized in two lattice sites (although 
there is a noticeable spread into neighboring sites which 
becomes more pronounced at lower U/t). This allows 
us to, for simplicity purposes, characterize the optimized 
solution in terms of a tiling scheme in the on-site basis 
(see bottom of the figure), although we emphasize that 
this is only approximate. In this case, the structure is 
defined by local dimers which are spin polarized (with a 
non-zero magnetization in each site) to yield an overall 
Neel-like configuration. 


B. ID: half-fllling 

We start by considering the half-filled ID periodic 
case. All calculations in this section, unless explic¬ 
itly stated, were performed in a periodic lattice with 
L = 120 sites, which we deem large enough to provide 
near-thermodynamic limit results for U/t > 1. Only uni¬ 
form tiling schemes were considered; clusters were defined 
in terms of a continuous set of I lattice sites, each filled 
with 1/2 electrons (for even 1). For U-cMF calculations 
with odd I, we have adopted a staggered configuration: 
if a cluster of size I has {I 1)/2 f-electrons and (Z — l)/2 
4 ,-electrons, its neighbors have {I -I- l)/2 4--electrons and 
{I — l)/2 t-electrons, respectively. As described in Sec. 
WM some spreading of the orbitals into neighboring 
sites is observed, particularly at low U/t. We note that 
broken-symmetry U-cMF solutions maintain the overall 
Neel-like structure observed in UHF, that is, a non-zero 


















guess orbitals optimized orbitals 



FIG. 1. (Color online) In the top schemes we show initial guess (left) and optimized (right) spin-orbitals that define a single 
cluster in an U-cMF calculation (12-site half-filled periodic ID lattice, Ujt = A) using clusters of 2 t- and 2 4.-orbitals with 1 
t- and 1 4,-electron. Orbitals are depicted (one per row) in the 12-site lattice (marked by small -|- signs) using the following 
conventions: t- (r-) orbitals are plotted in blue (red); filled (empty) circles indicate a positive (negative) orbital coefficient; the 
area enclosed by the circle is proportional to The guess orbitals (left) correspond to Boys-localized orbitals of the UHF 

solution (orbitals 1 and 2 are occupied; 3 and 4 are empty). The U-cMF optimized orbitals (right) displayed are those that 
diagonalize the t- and 4^- one-particle reduced density matrix, with orbitals 1 and 2 (3 and 4) having occupation of 0.9 (0.1). 
Note that the orbitals remain fairly local in character within a subset of 2 lattice sites. Accordingly, the bottom scheme shows 
a simplified representation of the optimized solution expressed in the on-site basis. This is characterized by local dimer-like 
structures, which we connect by a solid blue line. 


magnetization develops on each lattice site. 

We present in Fig. the energy per site obtained in 
cMF calculations at C//t = 2 (left) and U/t = A (right) 
as a function of the inverse of the cluster size, using both 
restricted (R-cMF) and unrestricted (U-cMF) optimized 
orbitals, as well as cMF calculations in the on-site ba¬ 
sis. We have also included, for comparison, the results of 
(exact) calculations carried out in a single cluster {L = I 
sites) using both open (OBC) and periodic boundary con¬ 
ditions (PBC). The former energies exactly match those 
of cMF calculations in the L — 120 lattice performed in 
the on-site basis, without orbital optimization. We note 
that L = I calculations using PBC do not provide a vari¬ 
ational estimate of the energy per site of the L = 120 
lattice (see Fig. [5] at U/t = 4, where the exact energy is 
approached from below). 

Comparing the results of cMF calculations that include 
orbital optimization with those in the on-site basis, it is 
evident that orbital optimization affords a significant im¬ 
provement in the variational estimate of the ground state 
energy. In addition, cMF calculations using unrestricted 
orbitals provide a sizable improvement over the corre¬ 
sponding restricted calculations at U/t = 4. We note 
that, at U/t = 4, cMF calculations do not converge (in 
l/l) as fast to the L — 120 limit as calculations using 
PBC, though the former have the advantage of being 
variational. On the other hand, at U/t = 2 cMF pro¬ 
vides better estimates than L = I calculations for t < 12, 
suggesting that a finite-size extrapolation with cMF re¬ 
sults should be preferred. We emphasize the significance 
of this given that, for arbitrary systems, an exact diago- 
nalization can currently only be performed up to lattices 
of size 18 or so. A linear extrapolation in l/l of U-cMF 
energies (using the I = 8, 10, and 12 results) yield the 
following estimates in the 1 —>■ L = 120 limit for the 
ground state energy E/{Lt): —0.5709(2) and —0.8414(2) 


for U/t = 4 and 2, respectively. These can be compared 
with the exact energies of —0.5738 and —0.8444. 

We will often refer to the fraction of correlation energy 
in assessing the quality of the ground state energy. The 
correlation energy is here defined as 

Acorr — Aexact AuHF: (^^) 

i.e., the difference between the exact and the UHF en¬ 
ergies. (Note that this differs from the traditional quan- 
tum chemistry definition based on restricted HF p^.) 
Figure [3] shows the fraction of correlation recovered in 
R-cMF and U-cMF calculations using clusters of increas¬ 
ing size {I = 2 to I = 12) as a function of U/t. The inset 
shows —Ecori/(Lt) as a function of U/t. The latter peaks 
at ~ 0.1 at U/t = 4. 

The fraction of correlation in restricted calculations us¬ 
ing I = 2 seems to vanish at large U/t, indicating that 
the Heisenberg limit predicted by this method is roughly 
the same as the UHF limit. This is not the case in un¬ 
restricted I = 2 calculations, which still recover around 
50% of the correlation in the large U/t limit. As I be¬ 
comes larger, the difference between restricted and un¬ 
restricted calculations gets smaller, as expected. U-cMF 
calculations with I = 12 are able to recover 85 to 90 % 
of the correlation energy across the entire U/t domain 
plotted, which spans both the weak and the strongly- 
correlated regimes. Accordingly, the maximum error in 
the energy per site in / = 12 U-cMF calculations is about 
0.011 at U/t = 4, which is remarkable given the simplic¬ 
ity of the approach. 

Figure 0] displays the fraction of correlation recovered 
in U-cMF and U-cPT2 calculations; results are compared 
with UMP2 and UCCSD. U-cPT2 energies significantly 
improve over U-cMF results for small cluster sizes. For 
instance, with I = 2 U-cPT2 recovers 90% (> 75%) of 
the correlation energy at small (large) U/t. Notice also 








9 





FIG. 2. (Color online) Energy per site obtained in cMF calculations on a L = 120 half-filled periodic ID lattice at Ujt = 2 
(left) and Ujt = ^ (right) as a function of the inverse of the cluster size 1. I = 1 R-cMF and U-cMF results correspond to 
restricted HF and UHF, respectively. We consider cMF results in the on-site basis and with a fully-optimized single-particle 
basis (opt). They are compared with (exact) calculations on a single cluster (i.e., L = 1) using open (OBC) and periodic (PBC) 
boundary conditions. The latter were computed by solving the corresponding Lieb-Wu [2^ equations. The energies oi L — I 
calculations with OBC exactly match those of cMF calculations in the full L — 120 lattice using the on-site basis. 
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FIG. 3. (Color online) Fraction of correlation (with respect to UHF) recovered in restricted (left) and unrestricted (right) cMF 
calculations in a L = 120 periodic ID lattice, as a function of U/t. (A log-2 scale is used in U/t for clarity purposes; results 
are shown from U/t = 1 to 16.) Cluster sizes from 2 to 12 were used. The inset in the right panel shows the total correlation 
energy per site, as a function of U/t. 


that U-cPT2 results do not overly deteriorate for large 
U/t as UMP2 does. UCCSD and U-cPT2 {I = 6) recover 
^ 90 % of the correlation energy across the entire U/t. 
It is interesting to point out that at C//t = 1 UMP2 and 
UCCSD results are better than U-cPT2 results with even 
Z; U-cPT2 results with odd Z, on the other hand, slightly 
overshoot the exact result. Unfortunately, the steep com¬ 
putational scaling of U-cPT2 rendered calculations with 
Z > 6 as too expensive with our current implementation. 

Figure [5] displays the spectrum of a single cluster 


Hamiltonian in U-cMF calculations, using cluster sizes of 
2,4, and 6. As PBC are used along with a uniform tiling 
scheme, all clusters end up displaying identical spectra, 
although we emphasize that this was not imposed. The 
eigenvalues are shown in Fig. [5] according to the and 
ni quantum numbers within the cluster, reference to the 
half-filled case. Here, the ground state in the (0,0) sec¬ 
tor of each cluster is used to construct the cMF state 
|$o). Although not shown in the figure, the spectrum 
should resemble, as Z becomes larger, that of a lattice of 
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FIG. 4. (Color online) Fraction of correlation (with respect to 
UHF) recovered in unrestricted cMF and cPT2 calculations 
m a. L = 120 periodic ID lattice as a function of U/t. UMP2 
and UCCSD results are provided for comparison purposes. 

I sites with OBC to the extent that orbitals remain fully 
localized. If denotes the ground state in the {x, y) 
sector, the perturbation series is stable (be., all denom¬ 
inators are positive) as long as the following conditions 
are met (we indicate in parenthesis the relevant interac¬ 
tions): 

1 . the (0,0) sector is gapped (two-cluster dispersion), 

2. 2eg°’°^ < (charge-transfer), 

3. 2eg°’°^ < (spin-flip), 

4. 2eg°’°^ < (2 clusters, two-electron 

charge-transfer). 

All the conditions are met in the cases plotted in Fig. [5l 
although as I becomes larger, 2eQ°’°^ ~ 

We show in Fig. [6] the contributions of different chan¬ 
nels to the second-order energy in U-cPT2 calculations as 
a function of Uft. Results from Fig. [5] indicate that two- 
cluster spin-flip (two neighboring clusters undergoing a 
spin flip) and two-cluster 1-electron charge transfer (two 
clusters interchanging a single electron) processes are the 
most important contributors to the second-order energy. 
Beyond that, the remaining two- and three-cluster inter¬ 
actions are small but non-negligible for U/t < 4. Four- 
cluster interactions are very small across all U/t. 

We now turn our attention to the spin-spin correla¬ 
tions in the ground state. As the cMF ansatz breaks the 
translational invariance of the Hubbard model, we have 
computed averaged spin-spin correlations, defined by 

'^(j) = ^E(Sj'-Sj-y), (36) 

j' 



{6n^,8n]) 

FIG. 5. (Color online) Spectrum of the cluster Hamilto¬ 
nian He in U-cMF calculations on a half-filled ID lattice 
at U/t = 4, using cluster sizes of 2, 4, and 6. All clusters 
in the L = 120 lattice display an identical spectrum. The 
eigenvalues are classified according to the and quan¬ 
tum numbers within the cluster. Here, (0,0) corresponds to 
a half-filled cluster, with 1/2 4- and ^electrons. Only those 
Hilbert sectors relevant to the evaluation of the second-order 
energy are shown. (Certain Hilbert sectors are missing as 
their spectrum is identical to one that is actually displayed; 
e.q., the (-1-1, 0) and (0, -1-1) spectra are equivalent, and so are 
the (-fl,-l) and (-1,-bl).) 
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loggF/t 

FIG. 6. (Color online) Contributions to the second-order en¬ 
ergy in U-cPT2 calculations (L = 120 ID periodic lattice) 
using a cluster size of 4 as a function of U/t. The notation 
used in the key takes the form clusters: interaction type”. 
Here, CT stands for charge transfer and SF refers to spin flip. 
The interaction channels are those from Tab. [H two-electron 
charge transfer processes are of opposite spin due to the na¬ 
ture of the Hubbard interaction. 


where j labels a lattice site. We plot in Fig. 0 the (real- 
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space) spin-spin correlations obtained from R-cMF and 
U-cMF calculations at U/t = 4. 

It becomes evident that unrestricted calculations (left 
panel) yield a structure with long range order. R-cMF, 
on the other hand, has non-vanishing spin-spin correla¬ 
tions only within the cluster; inter-cluster correlations 
vanish due to the spin singlet character of each cluster 
state. The long range spin-spin correlations in U-CMF 
are systematically decreased as the size of the cluster is 
increased. In the short range, both R-cMF and U-cMF 
yield significant corrections to RHF and UHF, respec¬ 
tively. With 1 = 8, R-cMF and U-cMF display simi¬ 
lar correlations to the first few neighbors (small k), with 
hints of the 1/k decay [H, present in the exact solu¬ 
tion at long range. 

A cleaner picture of the spin-spin correlations can be 
obtained by looking at them in reciprocal space. Figure[8] 
displays the (discrete) Fourier-transformed spin-spin cor¬ 
relations, at wave-vectors q = 0 and q = tt, obtained 
from R-cMF and U-cMF calculations. (R-cMF results 
at q = 0 are not shown as they identically vanish.) The 
two momenta are the most relevant ones: the q = 0 re¬ 
sult provides (S^), which should identically vanish for a 
spin singlet, while q = n provides the anti-ferromagnetic 
structure factor. We see that both values are significantly 
decreased in U-cMF with respect to what UHF predicts. 
Note that both R-cMF and U-cMF should converge to 
the same value (the exact one) as I tends to the L = 120 
limit. 


C. 2D: half-filling 

We now consider the periodic two-dimensional, half- 
filled square lattice. All the calculations in this section 
use a 12 X 12 periodic lattice, which should provide near¬ 
thermodynamic limit estimates for U/t > 2. In 2D, we 
are able to find a plethora of local minima (with respect 
to orbital rotations and coefficients in the cMF expan¬ 
sion) corresponding to different (approximate) tiling pat¬ 
terns. In principle, one could argue that, given a fixed 
number of clusters with their associated quantum num¬ 
bers (n-j-, nj,), an optimal solution (a global minimum) 
exists which minimizes the energy. We did not attempt 
to locate it but we did consider several uniform-like pat¬ 
terns that converge to different local minima. 

We show in Fig. the tiling patterns adopted in U- 
cMF calculations on the 12 x 12 lattice. A label used to 
identify each pattern is also provided in the figure. In 
most of the tilings displayed, a staggered configuration 
was chosen over an otherwise uniform tiling as it leads 
to lower variational energies. A staggered dimer config¬ 
uration is used in clusters of size 2. For clusters of size 
4, we discuss results with a square-based (staggered, 4s) 
tiling and a z-shaped (4z) tiling. We have considered a 
staggered configuration in terms of slabs (6si) and hats 
(6hi) in connection with clusters of size 6. Clusters of 
size 8 with a staggered slab (4 x 2, 8s) and a z-shaped 


(8z) configuration are used, which can be thought of as 
simple dimers of the considered size-4 clusters. We fi¬ 
nally also examined staggered plaquette configurations 
with clusters of size 9 and 12. 

We should point that, in constrast to Ref. 0, all the 
considered tiling patterns lead to a qualitatively correct 
description of the ground state character, i.e., all struc¬ 
tures lead to a ground state density with non-zero mag¬ 
netization in a Neel-type configuration. This is because 
the ground state of each cluster was independently opti¬ 
mized, without requiring that all the clusters share the 
same ground state. 

The left panel of Fig. [TUI shows the correlation energy 
(divided by the UHF energy) predicted by U-cMF as a 
function of U/t. (These can be contrasted with auxiliary- 
field quantum Monte Carlo (AFQMC) results shown in 
Fig. [TUI) Here, size-12 U-cMF predicts that the correla¬ 
tion energy is « 1.5% (~ 10%) of the UHF energy at 
U/t = 2 {U/t = 16). Large differences in the correlation 
energy predicted are observed as the clusters get larger, 
particularly at large U/t, which is unsurprising. A larger 
cluster, irrespective of its shape, tends to yield a larger 
correlation energy than a smaller one, but a few excep¬ 
tions are observed. When clusters of the same size are 
compared in different tiling schemes {e..g., 4s and 4z), 
we observe that the more compact the cluster is, the bet¬ 
ter the variational estimate for the ground state energy 
becomes at large U/t. Thus, the square tiling pattern in 
size-4 clusters yields a larger correlation energy than the 
z-shaped one. 

The right panel of Fig. |TU] shows the difference be¬ 
tween the double occupancy {D = ®ite 

predicted in U-cMF calculations with respect to that of 
UHF (which is shown in the inset). UHF overestimates 
(underestimates) the double occupancy at small (large) 
U/t. A relatively systematic improvement is observed in 
the double occupancy as the cluster becomes bigger. 

To show that other tiling patterns do not lead to funda¬ 
mentally different results, we show in Fig. [TT]the correla¬ 
tion energy predicted in U-cMF calculations with differ¬ 
ent tiling schemes using clusters of size 6. The additional 
tiling schemes (aside from those in Fig. O are shown in 
Fig. |TU] Several of them lead to approximately the same 
ground state energies. As previously discussed, the more 
compact the clusters are, the better the variational esti¬ 
mate of the ground state energy at large U/t. The same 
may not be true at small U/t. For instance, the lowest 
energies obtained at U/t = 2 corresponded to structure 
6z- 

We show in Fig. |TU]the ground state correlation ener¬ 
gies, divided over the UHF energy, obtained from U-cMF 
and U-cPT2 calculations as a function of U/t. Results 
are compared with UMP2, UCCSD, and AFQMC [^ . 
which can be deemed as numerically exact at half-filling. 
UMP2 displays the same behavior observed in ID, with 
the correlation energy vanishing for large U/t. It is ev¬ 
ident that the U-cMF results are not competitive with 
UCCSD or AFQMC, even with the larger clusters from 
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FIG. 7. (Color online) Real-space spin-spin correlations computed from unrestricted (left) and restricted (right) optimized cMF 
states for a L = 120 periodic ID lattice at U/t = 4. The I — 1 result in the left panel corresponds to UHF, while that in the 
right panel corresponds to restricted HF {i.e., a product of plane-wave states). 
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FIG. 8. (Color online) Fourier-transformed averaged spin- 
spin correlations {L = 120 ID periodic lattice at U/t — 4) in 
U-cMF and R-cMF calculations as a function of the inverse 
of the cluster size. Two q values are plotted, namely 0 and 
TT (note that the q = 0 curve is enhanced by a factor of 50 
for clarity purposes). U-cMF and R-cMF calculations should 
tend to the same q = tt finite value asl ^ L. The I —>■ L limit 
at q = 0 of U-cMF results should be 0 as the exact ground 
state corresponds to a true singlet state. 


Fig. [ini For instance, U-cMF using structure 6si recov¬ 
ers < 50% of the correlation energy across all U/t. On 
the other hand, U-cPT2 provides a sizable improvement 
over U-cMF results, with structure 6si capturing 85 % 
of the correlation energy atU/t = 12. This is not far from 
UCCSD, which recovers ~ 90 % of the correlation energy 
predicted by AFQMC at U/t = 12. (Finite size effects 


account for most of the difference between UCCSD and 
AFQMC at U/t = 2.) 

Figure [TT] displays the (discrete) Fourier-transformed 
averaged spin-spin correlations, at wave-vectors q = 
(0,0) and q = (tt, tt), obtained from U-cMF calculations 
at U/t = 8 in the 2D lattice. As it was done in the ID 
case, the spin-spin correlations are averaged as the wave- 
function ansatz breaks the translational invariance of the 
lattice. It becomes evident that as the cluster becomes 
bigger (regardless of the specific tiling pattern), the spin- 
spin correlations get reduced with respect to UHF. In 
particular, large clusters display less than half the spin- 
contamination per-site {i.e., the deviation of 5(0, 0) from 
0) of UHF. The anti-ferromagnetic structure factor is also 
reduced, though this should converge to a finite value in 
the limit I ^ L. 


D. 2D: lightly-doped regime 

In the 2D lightly-doped regime, we considered periodic 
square lattices with (n) = 0.8 and (n) = 0.875, following 
Ref. [13. We have used a 10 x 10 lattice for the former case 
and a 16 X 8 lattice for the latter case in order to have a 
lattice commensurate with the striped order expected to 
develop, as described below. 

Xu et al. [ 4 ^ discussed the UHF phase diagram for the 
2D Hubbard model in the half-filled and lightly-doped 
regime. For (n) ^ 0.9, the system transitions from a 
paramagnetic to a linear spin density wave regime at 
U/t 1 (c/. Fig. 17 in Ref. IdSh . A phase transition 
to a regime with diagonal spin density wave character 
occurs at U/t ss 4. Finally, the system becomes ferro¬ 
magnetic at large U/t. This UHF description has guided 
our U-cMF calculations. 
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FIG. 9. (Color online) Tiling patterns adopted in U-cMF calculations on a square 12 x 12 periodic 2D square lattice at half¬ 
filling. Sites within the same cluster are connected by solid lines; fully enclosed regions are shaded. A key is provided to the 
left of each structure. 


I. (n) = 0.8 

We show in Fig. [15] the spin and hole density profiles of 
UHF solutions with linear and diagonal spin density wave 
character. The (approximate) tiling pattern adopted in 
U-cMF calculations is superimposed in the figure. Clus¬ 
ters of size 6 with 4 electrons (2 f and 2 ),) have been 
used to describe the sectors with high hole density. On 
the other hand, the remaining regions with Neel char¬ 
acter have been described in terms of a staggered dimer 
configuration in structures 11 and Id, corresponding to 
the linear and diagonal spin wave character, respectively. 
In 21 and 2d, the dimers have been combined into half- 
filled clusters of size 4 as depicted in Fig. [151 


The resulting spin- and charge-density profiles ob¬ 
tained from U-cMF calculations with structures 1/ and 
Id are displayed in Fig. [151 The profiles resemble closely 
those obtained from UHF itself and shown in Fig. [T5] 
The main difference is the partial shift of the hole den¬ 
sity away from the main stripes into the neighboring sites. 
This suggests that the hole density is too localized in the 
UHF solution. 

We show in Tab. mi the resulting ground state ener¬ 
gies from U-cMF and U-cPT2 calculations. Results are 
compared with UHF, UMP2, UCCSD, and density ma¬ 
trix embedding (DMET) calculations [^, which can be 
deemed as highly accurate. U-cMF provides a significant 
improvement over UHF energies, indicating the inaccu- 
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FIG. 10. (Color online) (Left) Correlation energy (with respect to UHF), divided by the UHF energy, recovered in U-cMF 
calculations in a 12 x 12 periodic 2D square lattice, as a function of Ujt. The tiling patterns considered are those displayed in 
Fig- m (Right) Difference with respect to the double occupancy per site predicted by UHF in a variety of U-cMF calculations 
in a 12 X 12 periodic square lattice. The double occupancy of UHF itself is shown in the inset. 



logaU/f 


FIG. 11. (Color online) Same as Fig. HOI Different tiling pat¬ 
terns with clusters of size 6 (c/. Figs.[^and ll2D are displayed. 


rate treatment of short-range correlations in the simple 
HF description. The error in the UHF energies becomes 
very significant (> 0.21) at large U/t, sizably larger than 
the error in the half-filled regime. Interestingly, the lin¬ 
ear spin density wave character is favored in U-cMF even 
at relatively large U/t. We cannot rule out, neverthe¬ 
less, that this is an artifact of the particular tiling pat¬ 
tern chosen. The use of the size-4 clusters in place of 
the Neel dimers in U-cMF provides an improvement of 
about 10“^ t in the ground state energy for all U/t values 
quoted. 

Second-order PT provides a significant improvement 


over mean-field energies (both in HF and cMF). The 
UMP2 results are, however, far from UCCSD at large 
U/t. The difference between UCCSD(T) and UCCSD 
is also large in the strongly-correlated regime, indicat¬ 
ing the necessity of going beyond double excitations in 
the coupled-cluster ansatz. This is also evident by com¬ 
paring UCCSD and DMET results. U-cPT2(2/) is com¬ 
petitive with UCCSD at U/t = 4 but outperforms even 
UCCSD(T) in the large U/t regime. 


.2. (n) = 0.875 

We now turn our attention to even lighter doping, 
namely (n) = 0.875. We show in Fig. [T7]the spin and 
hole density profiles of UHF solutions with linear and di¬ 
agonal spin density wave character. Note that kinks are 
needed in a 16 x 8 lattice in the diagonal density wave 
profile; only a much larger 16 x 16 lattice is commen¬ 
surate with a fully diagonal profile. The (approximate) 
tiling pattern adopted in U-cMF calculations (with a lin¬ 
ear density wave character) is superimposed in the figure. 
As we did previously in the 10 x 10 lattice, the regions 
with Neel character have been described in terms of a 
staggered dimer configuration (structure 1 ), while the re¬ 
gions of high hole density are tiled into clusters of size 6 
(with 2 f- and 2 4,-electrons each). The Neel dimers have 
been combined into half-filled clusters of size 6 and 4 in 
structure 2, following the pattern indicated in Fig. [TTl 

Table nni shows the ground state energies obtained 
from U-cMF and U-cPT2 calculations. Results are com¬ 
pared with UHF, UMP2, UCCSD, and DMET. Just as 
in the (n) = 0.8 case, U-cMF improves significantly over 
UHF. It remains true that second-order PT provides a 
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FIG. 12. (Color online) Some additional tiling patterns (with clusters of size 6) adopted in U-cMF calculations on a square 
12 X 12 periodic 2D square lattice at half-filling. 


TABLE II. Ground state energies predicted with a variety of methods for a Hubbard 2D periodic 10 x 10 lattice with (n) = 0.8. 


method 

U = 2t 

II 

U = 6t 

U = 8t 

U = 12t 

DMET_^ 

-1.3062(4) 

-1.108(2) 

-0.977(4) 

-0.88(3) 


UHF (diag) 

-1.2165 

-0.9646 

-0.7933 

-0.6815 

-0.5501 

UHF (linear) 

-1.2678 

-0.9774 

-0.7843 

-0.6597 

b 

UMP22 

-1.3114 

-1.0760 

-0.8832 

-0.7767 

d 

UCGSD^ 

-1.3094 

-1.0925 

-0.9208 

-0.8246 

d 

uccsd7t)^ 

-1.3108 

-1.1045 

-0.9357 

-0.8444 

d 

U-cMF(ld) 

e 

e 

-0.8417 

-0.7396 

e 

U-cMF(2d) 

e 

e 

-0.8429 

-0.7406 

e 

U-cMF(lO 

e 

-T.0217 

-0.8520 

-0.7460 

-0.6271 

U-cMF(2Z) 

e 

-1.0227 

-0.8536 

-0.7478 

-0.6288 

U-cPT2(lO 


-1.0865 

-0.9380 

-0.8450 

-0.7394 

U-cPT2(2;) 


-1.0889 

-0.9435 

-0.8526 

-0.7513 


® Results extrapolated to the TDL from Refs.[^ and l43l. 

^ UHF fails to converge with this order. 

Calculations use the lowest energy UHF structure shown. 
Lower energy UHF solutions appear at large U/t. 

® U-cMF optimizations failed to converge. 


nice refinement on top of the mean-field result. The 
triples correction in UCCSD(T) becomes significant at 
large U/t, signaling the deficiencies in UCCSD. U-cPT2 
is a bit shy of UCCSD quality at U/t = A, but becomes 
competitive with UCCSD(T) at U/t = 8. 

Figures [IH] and [in] depict the Fourier-transformed spin- 
spin and density-density correlations obtained from U- 


cMF calculations (using structure 1). Just as in the 
case of ID, we have performed a global average over sites 
in order to remove the expected fluctuations due to the 
(spatial) symmetry broken character of the ansatz. The 
spin-spin correlations show a maximum at q = (Ttt/S, tt) 
which becomes more intense as U/t is increased from 4 
to 8. The density-density correlations display their max- 
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TABLE III. Ground state energies predicted with a variety of methods for a Hubbard 2D periodic 16 x 8 lattice with {n) = 0.875. 


method 

to 

II 

U = 6t 

U ^8t 

U = 12t 

DMET^ 

UHF (diag) 

-1.2721(6) 

b 

-1.031(3) 

b 

-0.86(1) 

-0.7184 

-0.6008 

-0.4682 

UHF (linear) 

-1.2270 

-0.9109 

-0.7128 

C 

C 

UMP2^ 

-1.2732 

-0.9858 

-0.7833 

-”6594 

e 

UCCSD^ 

-1.2719 

-1.0093 

-0.8305 

-0.7147 

e 

uccsdTt)^ 

-1.2738 

-1.0195 

-0.8446 

-0.7299 

e 

U-cMF(l) 

f 

-0.9476 

-0.7633 

-0.6477 

-”5176 

U-cMF(2) 

f 

-0.9500 

-0.7667 

-0.6514 

-0.5210 

U-cPT2(l) 


-1.0004 

-0.8289 

-0.7204 

-0.5965 

U-cPT2(2) 


-1.0040 

-0.8347 

-0.7276 

-0.6068 


^ Results extrapolated to the TDL from Refs. andl43l. 

^ UHF (diag) becomes UHF (linear) at low U/t. 

^ UHF fails to converge with this order. 

^ Calculations use the lowest energy UHF structure shown. 
® Lower energy UHF solutions appear at large U/t. 

^ U-cMF optimizations failed to converge. 



logatZ/i 


FIG. 13. (Color online) Same as Fig. 1101 U-cMF and U-cPT2 
calculations are compared with UMP2, UCCSD and AFQMC. 
AFQMC results, from Ref. [^ . correspond to thermodynamic 
limit estimates. 

imum at q = (7r/4,0). These observed profiles are con¬ 
sistent with the linear spin density wave character of the 
UHF charge and spin densities. Symmetry-projected cal¬ 
culations in Refs. 0 and [13 also show the same features. 
We refer the reader to Ref. [5l| for a discussion of the 
emergence of spin and charge order in the doped Hub¬ 
bard model. 


V. DISCUSSION 

In Sec. m we have described the cluster mean-field 
approach to treat strongly-correlated fermionic systems. 
A cMF state (including the orbital optimization degrees 
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FIG. 14. (Color online) Fourier-transformed averaged spin- 
spin correlations (12 x 12 periodic 2D square lattice a.t U/t = 
8) in U-cMF calculations as a function of the inverse of the 
cluster size. UHF {I = 1) results are also displayed. Two q 
values are plotted, namely (0,0) and (7r,7r). 


of freedom) is used as a variational ansatz for the ground 
state wavefunction. This, by construction, is guaranteed 
to provide better variational estimates than HF when the 
size of the cluster (assuming uniform tiling) is larger than 
1. Because of the simple nature of the ansatz, a RS-PT 
scheme can be adopted to account for the missing inter¬ 
cluster correlations. The results presented in Secs. IIVBI 
IIV Cl and IIV Dl provide evidence that a cluster-based ap¬ 
proach can provide a (semi)-quantitative description of 
the ground state of the half-filled ID and 2D Hubbard 
models, as well as for the lightly-doped regime in 2D 
square lattices. 


















17 


; fJ m ii 
: m ii 
i m li 
: fJ m ii 

J- -!■■■> -t-r^ > 'I' 


U 

a 

u 

u 

u 



FIG. 15. (Color online) Hole and spin density profiles of UHF solutions (10 x 10 periodic lattice, (n) = 0.8) obtained a.t U/t = 8 
with linear (left) and diagonal (right) spin density wave character. The magnitude of the hole density is proportional to the 
area of the green circles. The magnitude of the spin density is proportional to the size of the arrows. Superimposed on them 
we display the tiling patterns adopted in U-cMF calculations. Clusters with 6 orbitals and 4 electrons (solid blue) have been 
used in the high hole density sectors, while the Neel regions are described in terms of staggered dimers (solid red); this leads 
to structures 11 and Id, respectively, for linear and diagonal order. We have additionally considered a tiling pattern where the 
staggered dimers are combined into size 4-clusters, following the dotted lines, leading to structures 21 and 2d. 


•yy«#0yy0# 

•yy000yy00 

0yy0#0yy0# 

•yy0#0yy0# 

•yy0#0yy0# 

0yy0#0yy0# 

•yy0#0yy0# 

0yy0#0yy## 

0yy0#0yy0# 

• y ymm y y mm 


OO y y ooy 000 
QQ0 0 0 y ( 5 # y y 
y y##y yOQy 0 
0 y ooy 000^9 

000 yOO y y 0 O 
OO y y ® 0 y y y y 
OOy y y y## y y 
y y OO y 0 O 0 y 0 
00 9 % 0 00 yCO 
y y y 09900 99 


FIG. 16. (Color online) Hole and spin density profiles from U-cMF calculations (10 x 10 lattice, (n) = 0.8) at U/t = 8 with 
linear (left, 11) and diagonal (right. Id) ordering. The magnitude of the hole density is proportional to the area of the red 
circles. The magnitude of the spin density is proportional to the size of the arrows. 


In the half-filled ID model, results with comparable 
accuracy to UCCSD can be obtained by U-cMF using a 
sufficiently large cluster or by using U-cPT2 with smaller 
cluster sizes. Not only the energy is improved in U-cMF 
with respect to UHF, but also other ground state prop¬ 
erties such as spin-spin correlations. Due to the local na¬ 
ture of the interactions in the Hubbard Hamiltonian, con¬ 
tributions to the second-order energy arise mostly from 
two-cluster (spin flip and one-electron charge transfer) 
interactions. 

In the half-filled 2D square model, U-cMF was not as 
accurate as it was in the ID case, even when using clus¬ 
ters of size 12. This difference can be understood in terms 
of the missing inter-cluster correlations and the area-law 
of entanglement entropy [s^. Whereas in ID the size of 
the boundary (which determines the missing inter-cluster 
correlations) of a given cluster remains fixed, in 2D it 


scales as the perimeter of the cluster itself. This also ex¬ 
plains why more tightly packed clusters provide better 
energetic variational estimates for large U/t, as the clus¬ 
ters become more localized. A significant improvement to 
the ground state energy is obtained with U-cPT2, where 
results are again comparable (although slightly poorer) 
than UCCSD. 

In the lightly-doped regime, we have used ad hoc tiling 
schemes that mostly respect the underlying spin density 
wave profile obtained with UHF. Following this strat¬ 
egy, results will be relevant to the extent that UHF itself 
provides a qualitatively correct description of the char¬ 
acter of the ground state. Our calculations suggest that 
it does. U-cMF results provide a sizable improvement 
over the UHF description and U-cPT2 results are com¬ 
petitive with UCCSD at small U/t and with UCCSD (T) 
at large U/t. The predicted U-cPT2 energies are still 
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FIG. 17. (Color online) Hole and spin density profiles of UHF 
solutions (16 X 8 periodic lattice, (n) = 0.875) obtained at 
l//t = 8 with diagonal (top) and linear (bottom) spin den¬ 
sity wave character. The magnitude of the hole density is 
proportional to the area of the green circles. The magnitude 
of the spin density is proportional to the size of the arrows. 
Superimposed on the bottom scheme we display the tiling 
pattern adopted in U-cMF calculations. The Neel regions are 
described in terms of staggered dimers (solid red), leading to 
structure 1 , while regions of high hole density are described 
with slab-shaped clusters of size 6 (solid blue). We have also 
considered a pattern (structure 2) where the dimers are com¬ 
bined into half-Hlled clusters of size 6 and 4, following the 
dotted lines. 


above the DMET estimates, indicating that part of the 
long-range correlations expected to develop in the lightly- 
doped regime are still unaccounted for. These may be de¬ 
scribed with higher order RS-PT or other more powerful 
many-body approaches. 

We note that it is generally true that describing inter¬ 
cluster correlations (here via second-order RS-PT) im¬ 
proves the ground state energy to a larger extent than 
enlarging the size of the cluster in mean-field calcula¬ 
tions. (A clear exception to this is the fact that U-cMF 
with a cluster of size 2 provides better results at large 
U/t than UMP2.) The good quality of U-cPT2 results 
suggest that the renormalized Hamiltonian (expressed in 
terms of cluster states) is more amenable to a perturba¬ 
tive treatment than in the case of the HF particle-hole 
transformation. Thus traditional many-body approaches 
(such as second-order RS-PT) can be built on top of the 
cluster mean-field description to provide a high quality 
answer. This is further explored in Appendix]^ were we 
show the results of high-order RS-PT calculations in a 


FIG. 18. (Golor online) Fourier-transform S{qx,qy) of the av¬ 
eraged spin-spin correlations obtained from U-cMF (structure 
1) calculations in a 16 x 8 ((n) = 0.875) lattice at [//t = 4 
(top) and U/t = 8 (bottom). 
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FIG. 19. (Golor online) Fourier-transform M{qx,qy) of the 
averaged density-density correlations obtained from U-cMF 
(structure 1) calculations in a 16 x 8 ((n) = 0.875) lattice at 
U/t = A (top) and U/t = 8 (bottom). The strong peak at 
(0,0) (= /L) has been removed for clarity purposes. 
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small lattice. In the remainder of this section we discuss 
possible strategies that can improve the results presented 
in this manuscript. 

Perhaps the simplest strategy is to increase the flexi¬ 
bility in the mean-field variational ansatz, which can be 
done in a variety of ways. The full Hilbert space (be., not 
restricted to a given TOs sector) or, in fact, the full even- 
or odd-number parity Fock space within each cluster can 
be used. [ 5 ^ In doing this, it is not necessary to use a more 
general form for the single-particle transformation that 
dehnes the orbital optimization. The latter can be done 
in addition to (or in place of). Thus in systems where 
local number fluctuations are essential, a Bogoliubov-de 
Gennes single-particle transformation can be used. 

The local character of the clusters can be exploited 
in cPT2 calculations by, e.g., truncating the computed 
interactions according to some distance criterion. This 
could alleviate significantly the computational cost in 
cPT2 calculations (bringing them to linear scaling in the 
number of clusters) while also facilitating carrying out 
the cPTn expansion to a higher order. In order to deal 
with the large number of interacting clusters and states, 
a stochastic sampling of contributing processes can be 
performed. 

At this point, we would like to comment on the nature 
of the states used to carry out the perturbation expan¬ 
sion. In this work, we have used an energy criterion to 
truncate the number of cluster states when this was im¬ 
perative. Other criteria may be used, such as a density 
matrix based criterion (akin to the one used in DMRG 
[3Hi)- Here, one would diagonalize the Hamiltonian of 
the cluster interacting with part of its environment. The 
resulting ground state wavefunction is projected into the 
cluster states; those states with highest occupation con¬ 
stitute the optimal subset of states to use. Our main con¬ 
cern regarding this strategy is that the resulting cluster 
-|- (relevant) environment may become too large to solve 
(exactly) for its ground state. For instance, the environ¬ 
ment around a four-site square cluster in a 2D lattice 
should include at least eight additional sites/orbitals. 

Of course, other strategies to account for inter-cluster 
correlations may be used. One possible alternative in 
the case where there are a few nearly degenerate states 
in each cluster, is to diagonalize the full Hamiltonian in 
the direct product basis spanned by the relevant clus¬ 
ter states. This is part of the essence of the contrac¬ 
tor renormalization group (CORE) algorithm [53 - l3 | and 
has also been used in the active space decomposition 
(ASD) [3 [13 method in quantum chemistry. 

We think a coupled-cluster based Mproach such as 
the one proposed by Li in BCCC [3] is among the 
most promising avenues. In particular, a coupled-cluster 
ansatz should provide an improved description of the 
missing inter-cluster correlations in the mean field than 
low-order RS-PT. Given that the cluster-based Hamilto¬ 
nian contains up to four-tile interactions, it appears that 
the minimal coupled-cluster model should include up to 
quadruple excitations. Nevertheless, we have observed 


that two-tile interactions dominate the contribution to 
the second-order energy. It may not be unreasonable to 
restrict the excitation to singles and doubles. Moreover, 
locality can also be exploited within a coupled-cluster 
framework. 

Lastly, we would like to point out that even though we 
have used the cMF approach to study strongly interact¬ 
ing systems, it may be used in other contexts. In particu¬ 
lar, systems which can be effectively represented in terms 
of weakly interacting fragments of otherwise strongly- 
correlated fermions (see Appendix |B]) can be very effi¬ 
ciently described by low-order perturbation theory based 
on a cMF state. 


VI. CONCLUSIONS 

We have introduced a cluster mean-field variational 
approach and discussed its applicability to describe the 
ground state of strongly-correlated fermion systems. In 
this work, the full optimization of the cluster mean-field 
state has been carried out, including orbital optimiza¬ 
tion, with the restriction that the cluster state has well- 
defined n and rris quantum numbers. The restrictions are 
imposed in order to preserve N and Mg in the full sys¬ 
tem. The cluster product state constitutes an eigenstate 
of a mean-held (zero-th order) Hamiltonian, which allows 
us to formulate a RS perturbative approach to improve 
upon the mean-held description. 

We have presented mean-held and second-order per¬ 
turbative results of the ground state energies (and other 
observables) of the periodic ID and square 2D Hubbard 
models. In the half-hlled ID case, our U-cMF results 
become as accurate as UGGSD across all U/t for suffi¬ 
ciently large clusters. U-cPT2 results on smaller clusters 
also provide a consistent description across all interac¬ 
tion strengths. In 2D at half-hlling, U-cMF is poorer 
than in the ID case yet U-cPT2 provides ground state 
energies of near UCCSD quality. In the lightly-doped 
regime of the 2D model, U-cPT2 results remain competi¬ 
tive with UGGSD although they are still not competitive 
with DMET estimates. In general, we observe that U- 
cPT2 energies with small clusters are often better than 
U-cMF results with signihcantly larger ones. 

Overall, the results of this work suggest that a clus¬ 
ter mean-held approach can provide an excellent start¬ 
ing point and a path to a highly accurate, efficient de¬ 
scription of strongly-correlated fermionic systems, and 
the Hubbard model in particular. Several strategies to 
improve the mean-held description as well as correlated 
approaches built on top of it have been suggested. 
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Appendix A: High-order perturbation theory 

For sufficiently small systems, the standard RS per¬ 
turbation series can be evaluated to high order by direct 
solution of the RS-PT equations 

m —1 

1=0 

(Al) 

where we have assumed intermediate normalization 
((^f(o)|^M) = 0 V m > 0). We have evaluated the 
UMPn (ie., RS-PT using canonical UHF orbitals and 
orbital energies) and U-cPTn (as formulated in Sec. HIl 
using a cluster of size 2) perturbation series for a half- 
filled L = 8 ID periodic lattice at U/t = A and U/t = 8. 
The energy as a function of n is displayed in Fig. [20l 
The UMPn series approaches the exact energies very 
slowly, particularly at large U/t. On the other hand, 
U-cPTn is much faster approaching the exact energy, al¬ 
though the series has a divergent nature at t//t = 8. This 
is likely due to the near degeneracies expected to appear 
at large U/tin the spectrum of each cluster. It is possible 
that a convergent nature can be restored by tweaking the 
definition of the zero-th order Hamiltonian. In spite of 
that, these results support the premise that once correla¬ 
tions within the cluster have been described accurately, 
the ground state of the resulting renormalized Hamilto¬ 
nian can be expressed by a many-body expansion that is 
more rapidly convergent than the common UMPn series. 


Appendix B: Dimerized Hubbard model 

Consider a Hubbard lattice tiled into clusters. As the 
clusters become non-interacting, the cMF approach be¬ 
comes exact. In this section, we assess the quality of the 
cMF and cPT2 ground state energies as a function of 
the interaction between clusters. Being more specific, we 
consider the dimerized periodic Hubbard ID model (see, 
e.^., Ref. given by the Hamiltonian 

H = - h Y ^4+1, 

odd j,<7 

“^2 Y +h.c.)-f 

even j,(T j 

(Bl) 

Figure [m shows the UHF, U-cMF, and U-cPT2 ground 
state energies obtained in a half-filled L = 8 periodic 
ID lattice, as a function of the ratio t 2 /ti with U/ti = 
4. Clusters of size 2 have been used in cMF and cPT2 
calculations. 

At t 2 /ti = 0, the exact ground state energy reduces 
to four-times the energy of an L = 2 half-filled lat¬ 
tice with OBC. Naturally, U-cMF reproduces this result, 
while UHF converges to an energy that is equal to four- 
times the energy of the corresponding UHF result for 
the L = 2 lattice. U-cMF (U-cPT2) remains highly ac¬ 
curate up to t 2 /ti ~ 0.3 {t 2 /ti ~ 0.6). On the other 
hand, U-cMF and U-cPT2 are not nearly as accurate in 
the vicinity oi t 2 /ti = 1, yet they still provide sizable 
improvements over the UHF description. The mean-field 
approach recovers « 35 % of the correlation energy, while 
U-cPT2 is able to capture around 65 % of it. These re¬ 
sults suggest that a simple cluster-based approach can 
accurately describe weak interactions among otherwise 
strongly-correlated fragments. 
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